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Abstract 

We propose a novel approach for smoothing on surfaces, namely estimating a func¬ 
tion starting from noisy and discrete measurements. More precisely, we aim at esti¬ 
mating functions lying on a surface represented by NURBS, which are geometrical 
representations commonly used in industrial applications. The estimation is based 
on the minimization of a penalized least-square functional. The latter is equiva¬ 
lent to solve a 4th-order Partial Differential Equation (PDE). In this context, we 
use Isogeometric Analysis (IGA) for the numerical approximation of such surface 
PDE, leading to an IsoGeometric Smoothing (IGS) method for fitting data spatially 
distributed on a surface. Indeed, IGA facilitates encapsulating the exact geometri¬ 
cal representation of the surface in the analysis and also allows the use of at least 
globally —continuous NURBS basis functions for which the 4th-order PDE can 
be solved using the standard Galerkin method. We show the performance of the 
proposed IGS method by means of numerical simulations and we apply it to the 
estimation of the pressure coefficient, and associated aerodynamic force on a winglet 
of the SOAR space shuttle. 

Keywords: Functional Data Analysis; Isogeometric Analysis; Smoothing on Sur¬ 
faces. 

1 Introduction 

The estimation of a function from a set of noisy data is a very common task, which is often tackled 
by minimization of a penalized least-square functional, where the penalty involves a differential 
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operator, commonly based on second order derivatives, and occasionally on derivatives of a differ¬ 
ent order. Classical examples are offered by smoothing splines (see, e.g., Ramsay and Silverman 
(2005) and references therein) for estimating functions defined over real intervals, thin-plate splines 
(see, e.g., Wahba (1990), Wood (2006) and references therein). Soap film smoothing (Wood et al., 
2008), splines over triangulations (see, e.g., Lai and Schumaker (2007)) for estimating functions 
defined over regions of M^, and spherical splines (see, e.g., Wahba (1981); Baramidze et al. (2006); 
Alfeld et al. (1996)) for estimating functions defined over spheres or spheres-like surfaces. In this 
context, one of the main challenges in minimizing such penalized least-square functional consists 
in determining a suitable finite dimensional space representing, at the discrete level, the infinite 
dimensional space to which the function belongs. In other words, the challenge is to find a finite 
dimensional problem which is tractable and whose solution is close to the solution in the infinite 
dimensional space. The same challenge arises when looking for the numerical solution of a PDE. 
This common goal has been recently exploited to develop statistical tools to deal with spatially 
distributed data. In this respect, Ramsay (2002) considered planar smoothing optimizing a penal¬ 
ized least-square functional with a regularization term involving the Laplacian, and used the Finite 
Element Method (FEM) (see, e.g., Quarteroni and Valli (1994)) to solve the estimation problem. 
Sangalli et al. (2013) generalized the method proposed by Ramsay (2002) to include space-varying 
covariates and to account for boundary conditions, while Wilhelm (2013) explored a generalized 
linear version of the method. Azzimonti et al. (2014, 2015) further extended the model of San¬ 
galli et al. (2013) to account for any elliptic differential penalization, not necessarily based on the 
Laplacian operator. Ettinger et al. (2015) and Dassi et al. (2015) dealt with the case of functions 
defined over two dimensional Riemannian manifolds, by considering a regularizing term based on 
the Laplace-Beltrami operator and exploited a conformal parametrization of the manifold to solve 
an equivalent estimation problem on M^. Duchamp and Stuetzle (2003) also explored smoothing on 
complex surfaces using a penalization based on the Laplace-Beltrami operator. A common feature 
of all these contributions is the use of FEM to estimate the underlying function corresponding to 
the observed data. 

In this paper, we consider the estimation of functions defined over surfaces in starting from 
a discrete set of noisy observed data in points distributed on such surfaces. More precisely, we 
refer to sufficiently smooth functions defined over surfaces that can be represented by NURBS 
(Non-Uniform Rational Basis Splines). Indeed, NURBS are commonly used in Computer Aided 
Design (CAD) to represent most of the geometries of Engineering and industrial interest. From 
a more general point of view, the proposed model is a particular case of a Generalized Additive 
Model (GAM) (Hastie and Tibshirani, 1990) where the smooth component is defined on a surface. 
Hence, most of the theoretical results of GAMs can be directly applied to this model, e.g., for the 
quantification of uncertainty. 

IGA has been first introduced by Hughes et al. (2005) with the main idea of using the same 
basis functions to represent the geometry and then to approximate the solution of the PDFs defined 
in such computational domains. This facilitates encapsulating the exact geometrical representation 
when performing the analysis of PDFs defined in the computational domain; see Cottrell et al. 
(2009). In this respect, the isogeometric paradigm facilitates the use of an exact representation 
of the surface, while most of the current methodologies, as FEM, generally only handle its ap¬ 
proximation. This may induce an error on the solution due to the approximation of the geometry 
and may require complex meshing procedures. As mentioned, to estimate the function over the 
surface, starting from its discrete and noisy measurements, we minimize a least-square functional 
where the penalty involves the Laplace-Beltrami operator associated to the surface, analogously to 
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Figure 1: The SOAR shuttle highlighting the inboard winglet represented by a single NURBS 
patch. [Courtesy of S3, Swiss Space Systems Holding S'A], 


Ettinger et ah (2015), Duchamp and Stuetzle (2003) and Wahba (1981). The estimation problem 
is tackled directly on the surface and using IGA to numerically solve the associated surface PDE. 
Eor this reason, we name the resulting method IsoGeometric Smoothing (IGS). High-order PDEs 
can be straightforwardly solved by NURBS-based IGA; indeed, globally continuous NURBS 
basis functions can be easily defined on surfaces for some k = 0,...,p — 1, where p is the poly¬ 
nomial degree. This allows to use standard Galerkin method for numerically solving the PDE. In 
this respect, Dede and Quarteroni (2015) studied the use of IGA for 2nd-order PDEs defined on 
surfaces, Tagliabue et al. (2014) analysed IGA for high-order PDEs, while Bartezzaghi et al. (2015) 
for high-order PDEs defined on surfaces. All these works showed the efficiency and accuracy of 
IGA for the spatial approximation of surface PDEs. Partially related to our work, NURBS-based 
and IGA approaches are used in Beaubier et al. (2014) and in Dufour et al. (2015) as calibration 
tools. 

In this paper, we introduce the geometrical background to define the application framework 
of IGS. Then, we show how one can minimize a least-square functional involving a penalization 
term and hence estimate a function starting from a set of discrete and noisy measurements. We 
assess and compare IGS to the Thin Plate Splines (TPS) (Duchon, 1977; Wahba, 1990) by means 
of numerical simulations. Finally, as industrial application, we estimate the pressure coefficient 
and the aerodynamic force acting on the inboard winglet of the space shuttle SOAR, designed by 
S3, Swiss Space Systems Holding SA. S3 is a Swiss company currently developping, manufacturing, 
certifying, and operating a launch system for small satellites of weight inferior to 250 kg. A 
geometry of the SOAR suborbital shuttle for a preliminary study of aerodynamic forces is shown 
in Figure 1. In this application, the data are pointwise measurements of the pressure coefficient 
and the associated quantity of interest is the aerodynamic force. We propose a method to estimate 
functions defined over a NURBS surface starting from scattered noisy observations. 

This paper is organized as follows. In Section 2, we briefly recall B-splines and NURBS. In 
Section 3, we introduce some analytical tools and the IGS method. We show results for two 
numerical simulations in Section 4: the first one allows to compare IGS to TPS, while the second 
one corresponds to a complex NURBS surface, for which planar smoothers are not suited. Finally, 
in Section 5, we show the results of the estimation of the pressure coefficient and aerodynamic force 
over the winglet of the SOAR space shuttle. Gonclusions follow. 
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2 B-splines and NURBS 


NURBS are widely used in CAD for geometrical representation of surfaces (Piegl and Tiller, 1997). 
We start by defining an open knot vector of degree p as a sequence of values S = {.^i, ..., 
with knots < • • • < ^n+p+i, where the first and the last knots are repeated p + 1 times. The 
interior knots can be repeated at most p times and, if a knot is repeated m times, we say that its 
multiplicity is m. The B-spline basis functions Aj,p are defined recursively using the Cox-De Boor 
formula. We summarize some properties of B-spline basis functions of index i and degree p, 
for i = 1,..., n. 

• The basis function possesses p — m continuous derivatives across a knot of multiplicity 
1 < m < p and is C°°— continuous between the knots. 

• The support of the basis function Aj,p is compact and contained in p+1 knots spans [^j, . 

• The basis functions are pointwise nonnegative, i.e. > 0, for all i = 1,..., n. 

n 

• They form a partition of the unity, that is 5:W,p(0 = l,foralUE[ei,?„]. 

i=l 

Starting from the basis function Ni^p, the NURBS basis functions are defined as projective 

transformations of B-spline basis functions. Let wi,... ,Wn > 0 be positive weights, then, NURBS 
basis functions are defined as: 


RiAO = 


Wi 




NiAC), Vi = l,...,n. 


A NURBS curve C'(^) E with d = 2, 3 and control points Bi,..., E is defined as: 

n 

i=l 

To define surfaces with NURBS, we resort to the tensor product scheme. Given two open 
knot vectors S = {^i,... ,^„+p+i} and M’ = {ryi,... QiAO^ for i = 1,..., n, Mj^g{T]), 

for j = 1,... ,m, the corresponding univariate basis functions, the control points Bjj E M'^, and 
positive weights Wij for i = 1,..., n, j = 1,..., m, we define the bivariate NURBS basis functions 
as: 


Ri,jiC,r]) = 


QkAOMiAAwki'- 


and a NURBS surface as: 


1=1 j=i 

where D = (Cij ?m) x ivi^Vm)- For simplicity, we consider henceforth the same polynomial degree p 
along both the parametric directions, for which we rewrite the bivariate basis functions simply as 
Ri^piC, ■?/)) for i = 1,..., N^, where = nm. For an exhaustive description of NURBS, we refer 
the reader to Piegl and Tiller (1997). 


4 





0.0 0.2 0.4 0.6 0.8 1.0 



Figure 2: B-spline basis functions for different knots vectors. Basis functions of degree p = 2 and 
globally continuous (top) and of degree p = 3 and globally —continuous (bottom). 


We remember that in IGA, one uses the same basis functions to represent the surface and 
to approximate the solution of the PDF defined in such computational domain. In general, it 
is possible to enrich the NURBS basis without changing the geometry, i.e. by preserving the 
geometrical mapping, with the goal of obtaining a more accurate solution of the PDFs. In this 
respect, /i—refinement indicates a uniform knot insertion, which adds new basis functions while 
preserving the geometrical mapping, while p—refinement refers to an elevation of the polynomial 
degree of the basis functions, similarly to FFM. In addition, the so-called k— refinement is peculiar 
of NURBS and refers to a consecutive order elevation and a knot insertion which allow to increase 
the polynomial degree and continuity of basis functions. All the refinement procedures are discussed 
in details in Cottrell et al. (2009), while for an introduction to NURBS in the context of IGA, we 
refer the interested reader to Hughes et al. (2005). 


3 The IGS method 


We describe the mathematical framework of the IGS method. First, we introduce the surface 
differential operators and then we introduce the IGS method for smoothing functions on surfaces. 
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3.1 Geometrical framework 


Following Dede and Quarteroni (2015), let 12 C be an open and bounded parametric domain of 
finite measure with respect to the topology of Then, let S C be a compact, connected, and 
oriented surface, defined by a NURBS geometrical mapping X : 12 —>■ S such that: 

X : 12 C —)• S C s = (si, S 2 ) x = (xi, X 2 , X 3 ). (3.1) 

We assume that X is sufficiently smooth, e.g. X E (7^(12). Then, we define the Jacobian of the 
mapping X, denoted by VX as: 

dX 

VX:12^m3x2^ s^VX(s), (VX)ij(s) = ^(s), i = 1,2, 3, j = 1, 2. 

OSj 

We denote by (VW)i the i-th column of the matrix XX. The metric tensor of the mapping X is 
represented by: 

G : 12 ^ s ^ G(s), G(s) = (VX(s))^ VX(s) 

and we denote with ^(s) the square root of the determinant of the metric tensor of the mapping X: 

s^g{s), g{s) = -^det(G(s)). 

The metric tensor of the mapping X is assumed to be invertible almost everywhere in 12. 


3.2 Differential operators 

Let (j) E C'^(S) be a smooth function defined over the surface S, i.e.: 

x^(/.(x), (/.eC2(S), 

such that we can define the differential operators on the surface S. Then, the projection operator 
P(x) on the tangent plane to the surface at the point x E S is given by: 

P(x) = (I - ns <8) ns)(x) = I - ns(x)ns(x)'^, Vx E B, 

where ns(x) is the unit normal vector to the surface at the point x and I is the identity matrix^. 
Then, the surface gradient operator, denoted by Vs, is defined as the projection of the gradient of 
a function extended in a tubular region containing B, i.e.: 

Vs(?2(x) := P(x)V0(x) = V(/)(x) - (ns(x)^V(/>(x))ns(x). 


The Laplace-Beltrami operator, which is the surface Laplacian operator, is expressed as: 


As(?i(x) = Vs • (Vs0(x)) = trace[P(x) V^0P(x)], Vx E B, 


where Vs • v(x) = trace(Vsv(x)), for all v E C^(B), is the surface divergence and 

(V^(/))jj(x) = g^.g^ (x) is the Hessian matrix of </>. These operators can be rewritten in the 

parametric domain 12 using the geometrical mapping (3.1) as: 


Vs</>(x) =VX{s)G-\s)V{^oX){s) and As(/>(x) 


9{s. 


-V 


g(s)G-^(s)V(<^oA)(s) 


(3.2) 


where s = A ^(x). 


^We compute nsix) as nslx) 


[VX]i(b)x[VX] 2 (s) 

||[VX]i(s)x[VX]2(s)||2 


where s = A ^(x). 
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3.3 Mathematical model 


Let us consider N points pi, ... ,pjv located on S for which the observed values are yi,... ,yN 
respectively. We assume that: 


Vi = f{Pi) + £i, y i = l,...,N, (3.3) 

where / is a sufficiently smooth function defined on S that we aim at estimating and et are 
independent observational errors of zero mean and constant variance. 

Given a positive smoothing parameter A, we aim at minimizing the following parameter depen¬ 
dent functional: 


N 

- v{pi))'^ + X (Asu)^ dS = ||y - vatH^ + X dS, (3.4) 

JS dE 


where y = (yi,..., yN)'^ and v^v = (u(pi),..., v{p]\f))'^ is the vector of evaluations of the general 
function v at the points pi,...,pAr. This functional is the same considered by Ettinger et al. 
(2015) and by Duchamp and Stuetzle (2003) that uses a finite element representation, and by 
Wahba (1981), that considers functions on spheres and proposes spherical splines. The functional 
(3.4) involves the surface Laplace-Beltrami operator, which is, roughly speaking, a measure of 
the curvature of the function related to the surface. The use of the Laplace-Beltrami operator 
ensures that the regularization is invariant to rigid transformations of the coordinate system, which 
is desirable both for theoretical and practical reasons. If the observational errors are normally 
distributed, the functional (3.4) can be viewed as a negative rescaled Gaussian penalized log- 
likelihood. Hence, in such case, minimizing (3.4) is equivalent to maximizing a penalized log- 
likelihood. Then, for a given positive parameter A, the estimation problem is: 

find / e .T : Ja(/) < Ja(u), Vu G (3.5) 


where T" is a suitable functional space to be defined to ensure the well-posedness of the problem 
and / is the estimate of / in the functional space J-, which should lie at least in i.e. 

F C Indeed, since / G H^(E) C (^^(S), the evaluation of / at the points pi,...,p 7 v is 

well defined. Since the problem (3.5) will be later associated to a PDE, some essential boundary 
conditions (Brezis, 1999) should be specified in relation with the choice of F. In the simpler context 
of functions defined over planar domains, and with the penalizing terms involving linear second 
order elliptic operators, the well-posedness of problem (3.4) is extensively discussed in Azzimonti 
et al. (2014) under different kind of boundary conditions (Dirichlet, Neumann, or Robin (see, e.g., 
Brezis (1999)). In the following, we prove the well-posedness of the estimation problem (3.5) in 
the particular case of homogeneous boundary conditions using the Lax-Milgram theorem, following 
Ramsay (2002) and Sangalli et al. (2013). Let Hq(S) be defined as: 

iLg(S) := G : Vsu • n = 0 and u = 0 on , (3.6) 

where n denotes the outward directed unit vector normal to the boundary of S. Then, one 
can characterize the solution of the minimization problem (3.5) and ensure the existence and the 
uniqueness of the solution in the case where / is assumed to lie in J-" = H^iF). With this aim, we 
recall the Lax-Milgram Theorem (see, e.g., Quarteroni and Valli (1994)): 
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Theorem 3.1 (Lax-Milgram). Let J- he a Hilbert spaee, H x H ^ M. a continuous and 

coercive bilinear form and F : H ^ M. a linear and continuous functional. Then, there exists a 
unique solution of the following problem: 

find u £ F : G{v,u) = F{v), '^v G F. 

Moreover, if G(-, •) is symmetric, then u £ F is the unique minimizer in F of the functional 
J : —)■ M, defined as 

J{v) = G{v,v) -2F{v). 

Lemma 3.1. Let v,u £ F = ifg(S) he two functions and A > 0 a positive smoothing parameter. 
Let G{-, •) and F{-) be defined as: 

= (vatjUat) + A / (Asf)(Asn) dS, and F(u) = (y, vat) , (3.7) 

JT. 

where vat = (?;(pi),..., t'(pAf))^ CLnd uat = (m(pi), ..., u{pn))'^ for some distinct points pi,..., pat 
on S. Then, the bilinear form G{-, ■) is coercive, continuous and symmetric and F{-) is linear and 
continuous in 77q(S). 

Proof. First, we note that | • defined as |f|H 2 (s) = dS, is equivalent to the norm 

II ■ IIh 2 (s) in This implies that there exists a constant > 0 such that \v\h‘ 2 ce) > 

do,r 2 ||'?^||_H' 2 ( 5 ]), \/v £ ddg(S) (see, e.g., Quarteroni (2015)). Then, we have: 

G{v,v) = (vat,vat) + A^(Asu)^ dS = ||vAr||i + X\v\h2^y:) ^ ^Co,n\\v\\jj 2 (^^), Vv G idg(S), 

and so G(-, •) is coercive. We now show the continuity of G{-, •). Since dd^(S) C (7*^(2), there exists 
a constant Gi^q^n such that ||vAr ||2 < C'i^o^Ar||t’||H2(s)) We have: 


G{v,u) = (vAr,UAr) + Ay^(Asf)(Asri) < ||vAr|| 2 ||uAr ||2 + A|z;|jA 2 (S)l'w|i^ 2 (S) 

^ C'f^O,7vll^ll//2(s)||^i||H2(s) + A|t>|jy2(s)|u|j^2(s) 

< max{Cl^^^,X}\\v\\H2(E')\\u\\H2(^j^), Vv, u G id|(S). 

Then, the bilinear form G{-, ■) is also continuous and its symmetry is obvious. 

Finally, we have: 

l-^(^^)l = I (yWAr) I < IlylbllvArIb < C'l,n,Ar||y||2||t’||H2(s), £ idg(S), 

which proves the continuity of the linear form F and concludes the proof. □ 

Proposition 3.1. Let F = ddg(S) and A > 0. Then, the solution of problem (3.5) exists and is 
unique. Moreover, problem (3.5) is equivalent to: 

find f £ F : ^vat, Jat) + A ^ A^u As/dS = (y, vat) , £ F. 


(3.8) 


Proof. The functional Jx{v) (3.4) can be rewritten as: 

J\{v) = ||y - VAr ||2 + A^(As/)^ dS = ||y ||2 - 2 (y, vat) + ||vAr ||2 + A^(Asu)^ dS. 
By defining Jx{v) as 

Jx{v) = llvArlla + A /" (Asu)^ dS - 2 (y, vat) , 
ds 

we have that 


arg min Jxiv) = arg min Jx{v). 


From the definitions of G{v,u) and F(v) of (3.7), the functional Jxiv) can be written as: 

Jxiv) = Giv,v) - 2Fiv). 


Thanks to the Lax-Milgram Theorem, it is then sufficient to use Lemma 3.1 to establish the well- 
posedness of problem (3.5). Moreover, since G'(-,-) is symmetric, problem (3.5) is equivalent to 
problem (3.8). □ 

Setting a priori the essential boundary conditions Vs/ • n = 0 and / = 0 on (9S following (3.6) 
may be an inadequate choice in several applications, especially when the behaviour of the function 
/ at the boundaries is not known a priori. In such cases, it may be more convenient to consider 
instead natural boundary conditions, that is: 


I Vs(As /) • n = 0 
[ As / = 0 

In the case that the boundary conditions (3.9) are set, we are unable to show the well-posedness of 
problem (3.8) with F = id^(S). Nevertheless, numerical experience indicates that it still yields a 
numerically stable problem. We can also note that the usual planar smoothers, such as TPS, also 
implicitly use natural boundary conditions. 

3.4 Numerical approximation: IGS 

Let be a finite dimensional approximation of / obtained by means of IGS. Let {'i/’i,..., 'f’xfh} be 
a basis of a discrete function space F F H’^iF) of dimension N^. In the finite dimensional 

space F^, problem (3.8) reads: 

find E F^ : (v^, {^) + X [ A^v^ As/'* dS = (y, , Vu'* E F\ (3.10) 

J s 

where v(y := (u^(pi),..., u^(pAr))^ and := (/'*(pi),...,/'*(p 7 v))^. Since F^ is finite dimen¬ 
sional, problem (3.10) is equivalent to: 

find/'*: <^i/)iAr, f/r) +^ As/'* As/'* dS = (y, i/jiAr) , Vi = 1..., M'*, (3.11) 


on 

on 


(3.9) 
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where := ('0i(pi),..., V’i(Piv))^- Let us define the x matrix R as 

(R)ij = AsV’* AsV'i ^ ^ matrix 'J' as (’J')ij = V'i(Pi)- Since belongs to 

it can be written as a linear combination of the basis functions: 

Afh 

= 'l2fi V’i(x), Vx e s, 

i=l 

or compactly as /^(x) = 'i/;^(x)f^, where := (fi,..., and 

■j/>(x) = ('01 (x),..., 0iY/i(x))'^. Then, we have: 

y:={^^ = 


Problem (3.11) in matrix form reads as: 

findf'^EM^^ Af^ = T'^y, 

where A := (T'^T' + AR). Then, the explicit form of is given by: 

f'^ = A'^T'^y = (T'^T'+ AR)“^T'^y (3.12) 

We see from (3.12) that the estimator has the typical form of a penalized least-square estimator. 
Since y = we finally get the evaluation of the function in the points {pi,..., pat} as: 

y = ^r(^r'r^r + AR)~^T''^y. (3.13) 

The smoothing matrix, that maps the observed data values y in the fitted data values y, is given 
by: 

S = T'(T'^T'AR)“^T'^. (3.14) 

The trace of S is a measure of the equivalent degrees of freedom of the estimator (see Buja et al. 
(1989)). If A = 0, the number of degrees of freedom is equivalent to the number of basis functions 
N^. However, the two notions differ for A > 0. While different definitions of equivalent degrees of 
freedom can be considered, these can be assumed as a consistent measure of the number of degrees 
of freedom that takes into account the harmonic penalization. In this respect, if the smoothing 
parameter A is strictly positive, the number of equivalent degrees of freedom is smaller than the 
number of basis functions used in 

We use IGA to solve the minimization problem (3.5), for which we define: 

= span o r?), i = 1,..., N’^} , 

where Ri^p are the NURBS basis functions used to build S, eventually after the application of 
some /i—, p— or /c—refinement procedure, as described in Section 2. is the number of basis 
functions, which is the dimension of J-^. NURBS allow to define basis functions which are globally 
C^—continuous on S. As consequence, one can approximate problem (3.5) with the standard 
Galerkin method, since C see Bartezzaghi et al. (2015) and Tagliabue et al. (2014). In 

this manner, we obtain a method, which we name IGS, allowing to perform smoothing on surfaces 
by means of NURBS-based IGA. This also allows encapsulating the original description of the 
geometry of the surface in the analysis. 
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The smoothing parameter A may be chosen by minimization of a generalized cross-validation 
criterion (GCV), defined as: 


GCV(A) 


N 

[N — trace (S(A))]^ 


l|y(A)-yf 


(3.15) 


see Graven and Wahba (1978). Here, we use trace (S(A)) as a measure of the equivalent degrees 
of freedom (EDF) of the model (Buja et ah, 1989). In order to solve the optimization problem 
corresponding to GGV minimization, we use a BFGS quasi-Newton method (see, e.g., Nocedal 
and Wright (1999)) with a sufficiently small tolerance. One can observe that the computation 
of the GGV criterion involves the inversion of the matrix + AR of size x N^. In our 

implementation we use a direct method to compute the matrix S because of the moderate size of 
this matrix. Methods based on matrix decompositions can improve both the efficiency and the 
stability of the optimization procedure used for the GGV criterion, (see, e.g.. Wood, 2006, pp. 
178-181). Other methods for the choice of the smoothing parameter are also available, e.g. the 
restricted maximum likelihood estimation (Wood, 2011). 

We remark that our model only considers a deterministic location of the measurement points, 
according to (3.3). While it is easy to account for random location of points at the implementation 
level of the IGS method, the uncertainty quantification in this setting is not straightforward. 


3.5 Distributional properties and quantification of uncertainty 

Here, we denote by E[y] and var(y) the expectation and the variance of y respectively. Moreover, 
let cr^ be the constant variance of the noise introduced in (3.3). For a given smoothing parameter 
A, the estimate is a linear transformation of the observations y, as shown in (3.12). Moreover, 
we have E[y] = i]\f = (/(pi), ..., f{pn))'^ and var(y) = cr^I. Then, from (3.12), we get: 

E[f'^] = E[(T'^T' + AR)^^T'^y] = 

Similarly, we can directly express the variance as: 

var(f'^) = var(A“^T''^y) = o-2A“^T''^T'A”\ (3.16) 

where we used the fact that the matrix A is symmetric and hence A~^ = A“^. In particular, 
under the assumption of normality of the errors, we have: 

Then, we can also express the expectation and the variance of the fitted values y explicitly as: 


E[y] = S fTv, 


and: 

var(y) = cr^SS^ = S^, 

respectively, since the smoothing matrix S is symmetric. In practice, the error variance must 
be estimated from the data. Following Hastie and Tibshirani (1990) we can estimate by: 


l|y -y|| 

N — trace(S) 


(3.17) 
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Given an additional point pat+i on S, the predicted value of the function / is given by: 


f^{PN+l) ='^fi MPN+i) = ^A(pAr+l)'^f^, 
i=l 

while its variance is given by 

var(/^(pAr+i)) = ^(pAr+i)'^var(f'*)r/)(pAr+i) = V(PiV+i)- 

An estimate of var(/^(pAr+i)), say vaf(/^(pjv+i), reads: 

^(/^(PAT+i) = (T^^(pAr+i)^A^^'J'^'J'A~V(PAf+i)- 

These results fully characterize the estimates in the case of Gaussian noise. In such case, one 
can derive confidence bands on the estimated functions and thus quantifying the uncertainty of the 
estimations of any predicted value. If the Gaussian assumption does not hold, this confidence inter¬ 
val should be used with caution and the confidence level is only approximated. More generally, the 
quantification of uncertainty has been widely studied in the context of generalized additive models 
(see, e.g.. Wood (2006) and references therein). A Bayesian approach to uncertainty quantification 
for this class of models is also possible; see Marra and Wood (2012). 

4 Numerical simulations 

In order to assess the IGS methodology, we consider two simulations on surfaces represented by 
NURBS for which the function / to be estimated is given a priori. We aim at showing different 
properties of IGS in different settings. In the first simulation, the configuration is such that any 
method used for planar smoothing such as Thin Plate Splines (TPS) (Duchon, 1977; Wahba, 1990) 
can be efficiently used and thus compared to IGS. In the second simulation, the setting is such that 
methods for two dimensional smoothing are not appropriate, while IGS can be straightforwardly 
used. 

4.1 Simulation 1 

We consider a quarter of cylinder S defined as S = = 1, 0 < z < 2, x > 0, y > 0}. 

Using cylindrical coordinates, it is parametrized by the following mapping: 

A : U = (0,2) X ^0, 0 S, A(si, S 2 ) = (cos(si), sin(si), S 2 ). 

This is an isometric mapping (Stoker, 1989). That means that the mapping preserves lengths and 
angles, and thus the area. In other words, the parametric domain is not distorted by the mapping. 
This kind of mapping allows to indifferently work on the parametric domain H or directly on the 
surface S. Thus, the smoothing can be performed on the parametric domain, which is planar, 
namely using any traditional method for planar smoothing. In particular, in the following, we shall 
compare the proposed technique to TPS. 

The surface S is exactly representable using NURBS basis functions of degree 2 or higher which 
are at least globally continuous. We are interested in recovering the function: 



f{x,y,z) = sin 
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Figure 3: Simulation 1. Surface S (quarter of cylinder) and exact function /. 


from noisy and discrete observations. The quarter of cylinder S and the function / are shown in 
Figure 3. The function / is evaluated in = 100 points located on a 10x10 grid of equally 

spaced points in the parametric domain and is affected by independent Gaussian observational 
errors £i of zero mean and standard deviation a = 0.125. We generate the data Ui = /(pi) + Cj 
for i = 1,. .., iV and the simulation is repeated M = 100 times. We use random scalars drawn 
from the standard normal distribution to generate the noise (specifically, we used the MATLAB 
function randn). The dimension of the IGA space varies between = 49 and = 121 
and we use different NURBS basis functions, namely globally C^ — , , and C^— continuous of 

degrees p = 2,3, and 4 respectively, obtained from /c—refinement. The smoothing parameter A 
is chosen at each simulation repetition and for each basis setting considered, by minimizing GGV 
criterion in (3.15). 

The first two estimated functions over the M ones are displayed in Figure 4. We observe 
that these are qualitatively good estimates of /. In general, the number of basis functions must 
be chosen carefully. Indeed, this choice should depend both on the complexity of the function to 
be estimated / and of the number of data points N available. When the number of basis functions 

is small, IGS is not able to capture the behaviour of the function /, as it would be the case 
with any other smoother. On the contrary, when the number of basis functions is high, we see 
that there is a larger variability in the estimated functions /^. Indeed, when the number of basis 
used is too high, GGV criterion can lead to overfitting, that is the estimated function incorporates 
noise. However, we see in Figure 4 that the estimates are not very sensitive to this choice. Finally, 
we remark that the minimum number of basis functions is basically dictated by the number of 
functions used to represent the surface with NURBS. 

We notice, following Section 3.3, that we used natural boundary conditions, for which we have 
not formally proved the well-posedness of the problem. We report in Table 1 the mean condition 
number for the matrix A of (3.12), with A chosen with the GGV criterion and for p = 2 only, 
since results for p = 3 and p = 4 are similar. As a matter of fact, the system (3.12) results to be 
well-conditioned in all our numerical experiments. In order to better assess the IGS method, we 
use the empirical mean function ^ associated empirical variance function 

^ “ /^(x))^, where denotes the i-th estimated function. These are both shown in 





(a) = 49 




V. 


(c) iV'* = 64 


(d) = 64 



Figure 4: Simulation 1. Estimated functions in the first (left) and second (right) simulation 
repetition, out of M = 100 repetitions, using = 49 (top), = 64 (middle) and = 121 
(bottom) number of basis functions, globally C^—continuous of degree p = 2. The corresponding 
measured values {yi}^i are displayed on the same color scale as the true function and estimates. 




V. 
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Figure 5: Simulation 1. Empirical mean function (left) and empirical variance of the function 
(right) over the M=100 simulation repetitions, for NURBS basis functions of degree p = 2, 
globally C^— continuous and of dimension = 81. 



25 

36 

49 

64 

81 

100 

121 

144 

169 

196 

Roo(A) 

46.1 

47.3 

48.2 

46.5 

43.9 

80.9 

210 

431 

355 

307 


Table 1: Simulation 1. Mean condition number R'oo(A) for the matrix A, with smoothing parameter 
chosen with GCV, p = 2 and for different number of basis functions N^. 

Figure 5, in the case of globally continuous basis functions of degree p = 2 and with = 81 
basis functions. The estimates appears to have a negligible bias and a relatively small variance. 

We compare our methodology with a widely used smoothing technique, TPS, using cylindrical 
coordinates. This method is implemented for instance in the R package mgcv (Wood, 2015). We 
use different number of basis functions for a comparison, selecting the smoothing parameter at 
each simulation repetition and for each number of basis considered via GCV. As a criterion for the 
comparison, we use the mean squared error (MSE) of the estimator, computed as: 

i=i 

where li,... ,1/, is a lattice of 150 x 200 evaluation points on S. We compute the MSE for = 
49,64, 81,100, and 121 and degree p = 2,3, and 4 for IGS and for = 40, 50,..., 100 for TPS. 
The comparisons of MSE in terms of the equivalent degree of freedom is shown in Figure 6. We also 
compare in Figure 7 the best setting of each methodology and for each number of basis functions, 
that is the lowest median MSE. Finally, Figures 6 and 7 show that IGS is comparable to TPS in 
terms of performance, in a setting where the latter technique may be applied. Figure 6 illustrates a 
key feature of the IGS method and more generally of all smoothing techniques. Specifically, one can 
notice that increasing the number of basis functions does not improve the estimated function /^. 
Indeed, although we increase the number of basis functions to build /^, the number of measurement 
points N remains the same, i.e. is still built from the same set of data, but only through a 
richer finite dimensional space J-^. Therefore, the convergence of to / should be simultaneously 
regarded through the number of data N and the quality of the NURBS space We highlight 
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Figure 6: Simulation 1. Comparisons of the MSE for IGS and TPS (top-left), in terms of the equiv¬ 
alent degrees of freedom (EDF). IGS uses globally C^—, and continuous basis functions 
of degrees p = 2 (top-right), 3 (bottom-left) and 4 (bottom-right), respectively. 


in Figure 6 that the same happens with TPS-based smoothing. A key role in the convergence of 
the method may be played by the NURBS space J-^ when the number of measurement points is 
clustered in a region of S, for which mesh refinement techniques can be used. 
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Figure 7: Simulation 1. Comparisons of the MSE for the best setting of IGS and TPS. 


4.2 Simulation 2 

We consider the surface S reported in Figure 8 which is represented in terms of NURBS basis func¬ 
tions of degree p = 2 starting from the knot vector S = {0,0,0,1,1,1} along both the parametric di¬ 
rections and control points Pi = (1,0,0)^, P2 = (0,0.75, 0)^, P3 = (0,1, 0)^, P4 = (1, 0,1)^, P5 = 
(1,1,0.S)"^, Pe = (0,1,1)^, P7 = (0.25,0, l)"^, Ps = (0.25,0.25,0.25)'^, and Pg = (0.25, 0,1)"^, and 
with the corresponding weights vector w = (1, l/\/2,1, l/\/2,1/2, l/\/2,1,1/v^, 1)^. As reported 
in Figure 8, we consider the exact function to be estimated: 

/5 

f{x, y, z) = (2 (xsin(7ry) -|- cos(7rx)y) — 1) cos ( -ttz 

which is evaluated in N = 100 points {pij^i located on S; more specifically, they are located on 
a 10x10 grid of equally spaced points in the parametric domain. The N observations are affected 
by independent Gaussian observational errors of zero mean and standard deviation a = 0.125. The 
sampling of the data y* = /(pi) -|- e*, i = is repeated M = 50 times. The sampling 

locations {paili remain the same for all the repetitions. We consider three NURBS spaces 
of dimensions = 100, of globally C^—, C^— and C^— continuous NURBS basis functions of 
degrees p = 2,3, and 4, respectively. These are obtained by fc—refinement of the original NURBS 
basis. 

In Figure 9 we report the estimated functions for the last two simulation repetitions, by 
considering NURBS spaces of basis functions of degree p = 3 and globally C^—continuous on 
S. Figure 10 highlights the empirical mean function and the corresponding empirical variance 
function over the M = 50 simulation repetitions, obtained for the same NURBS basis functions. 
The results obtained for the degrees p = 2 and 4 are very similar and are not reported here for the 
sake of brevity. Our estimation has a negligible bias and a small variance, as shown in Figure 10. 
The behaviour of the function / seems to be very well captured by our estimator /^, as shown in 
Figure 9. In this setting, we increase the regularity of the basis functions without changing the 
number of basis functions N^. The estimated functions are thus globally continuous. We 

then compare the mean integrated squared errors (MISE) of IGS for different regularity and degrees 
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Figure 8: Simulation 2. Surface S and exact function / 




Figure 9: Simulation 2. Functions estimated for the last two repetitions, for = 100 points; 
positions of the points and the corresponding measured values {yi}^^i for NURBS basis 

functions of degrees p = 3 and globally continuous on S. 


of NURBS basis functions. The MISE of any estimated function / is defined as: 

MISE(/) = ^/^ (/-/)' dS. 

As we observe in Eigure 11, the quality of the estimation is not significantly affected by the regularity 
of the basis functions. 
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Figure 10: Simulation 2. Empirical mean function (left) and empirical variance of the function 
jh for N = 100 and M = 50 obtained with NURBS basis functions of degrees p = 3 and 

globally continuous on S. 
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Figure 11: Simulation 2. Boxplots of the MISE for IGS with degrees p = 2,3, and 4, globally 
— and (7^—continuous basis functions. 


5 Estimation of aerodynamic force on the SOAR’s winglet 


We aim at estimating the pressure coefficient field Cp and the corresponding aerodynamic force 
on the inboard winglet S of the SOAR shuttle shown in Figure 1. The pressure coefficient is a 
dimensionless field related to the pressure field. It describes the relative pressure over the surface 
of the winglet and is given by: 


^ _ Pi^) - Poc 

- Poo „2 ■ 
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Figure 12: Original data (left) and smoothed estimated pressure coefficient C^, with N = 824 
data points and = 584 basis functions (right). 


where p(x) is the pressure at point x, Voc and poc are the far field wind speed and pressure, with 
Poo the air density. The flow regime of the free stream is subsonic, namely the Mach number is 
0.7, for which we can reasonably assume that the pressure coefficient field over the winglet remains 
sufficiently smooth, since transonic effects are marginal. The force F acting on the winglet and due 
to the contribution of the pressure (see, e.g., Anderson (2010)) is given by: 

^ = J^(P~ ^ 

where ns is the unit normal vector to the surface S. The final quantity of interest of the estimation 
is the aerodynamic force F. The pressure coefficient Cp is measured in = 824 data points on 
the surface S for which the sampled data are depicted in Figure 12, which shows the data after 
application of an affine transformation^, that modifies the values of about 3% at most. These data 
are derived from Computational Fluid Dynamics (CFD) simulations and represent a preliminary 
study prior experimental campaign in a wind tunnel. The inboard winglet is represented by a 
single NURBS patch surface built with = 584 degrees of freedom. The geometry is built with a 
NURBS basis of degree p = 5 and functions globally continuous. Thus, the minimum number 
of basis functions representing the function space for IGS is = 584. 

First, we estimate the pressure coefficient field Cp with all the available points. We use IGS to 
estimate the pressure coefficient field and then the aerodynamic force. By using IGS with N = 824 
and = 584, for NURBS basis functions of degree 5 and globally continuous basis functions, 
we obtain the Cp field reported in Figure 12 (right). 

We now assess the quality of IGS under the hypothesis that less data points than those effectively 
available can be used. Indeed, in industrial applications, the experimental measurements can be 
quite expensive and it could be impractical to have a large number of sampling points (in this case, 
pressure probes), and the pressure field can be measurable in fewer points than with numerical 
simulations. To assess the validity of IGS when few data points are available, we compute the 
aerodynamic force estimated from a subset of points over the winglet. We compare the results 

^The displayed data have been modified for copyright reasons with respect to the original ones provided by S3, 
Swiss Space Systems Holding SA. 
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N 

difference of direction (in degrees) 

||F'^-F((,f||/||F(*,f||(in%) 

584 

103 

1.007 

3.293 

584 

206 

0.816 

1.036 

584 

412 

0.084 

0.745 

584 

618 

0.001 

0.003 

989 

412 

0.278 

0.864 

989 

618 

0.061 

0.205 

989 

824 

0.003 

0.117 

2079 

824 

0.0113 

0.296 


Table 2: Relative difFerence between vector of forces with respect to the reference force F^ef) 
estimated in the case where = 584 and N = 824, in degrees (left) and in norm (right). Direction 
(in degrees) and relative modulus change (in %) of the estimated aerodynamic force obtained 
with different points N and basis functions with respect to the reference force F(*g£. 


in terms of the aerodynamic force, which is the quantity of interest, by using fewer data points and 
different number of basis functions with respect to the reference setting corresponding to = 584 
and N = 824. We compare the forces in terms of the relative difference of the modules (in %), that 
is and of the direction (angle in degrees). The results are reported in Table 2. 

We first use all the N = 824 available points depicted in the Figure 12 (left) to estimate the pressure 
coefficient field. The estimation of the reference force is given by F(ig£ = (0.5080,49.40,-26.03)^ 
kN, which is in line with the expectations in terms of direction and magnitude, according to the 
CFD simulations. 

We then use only N = 618,412,206, and 103 points and NURBS spaces of dimensions = 
584, 989, and 2079. The data sets with N = 103, 206 and 412 have been built using one, two and 
four over eight of the points in the original sequence, respectively; for N = 618, the points are the 
complementary of those in the set with N = 206. In Figure 13, we report the estimated pressure 
coefficient fields on the surface computed using IGS. The choice of the smoothing parameter is done 
via minimization of the GCV criterion. The estimated functions are very similar when the number 
of points is relatively large {N = 824, 618, and 412), as shown in Figure 13. The relative change in 
the estimated force is quite small, as we can observe on Table 2. IGS is able to accurately estimate 
the force even with a small subset of the original data. 
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-0.95 0.26 -0.95 0.26 


(c) N = 412 and = 594. 


(d) N = 618 and iV'* = 594. 



-0.95 0.26 -0.95 0.26 


(e) N = 824 and = 989. (f) N = 824 and = 2079. 

Figure 13: Estimated pressure coefficient field Cp based on different number of data points N and 
number of basis functions N^. 
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6 Conclusions and discussion 


In this paper, we proposed a methodology, IGS, to deal with functional data defined on surfaces 
described by NURBS, specifically to reconstruct the functions from noisy observations. The pro¬ 
posed IGS method has the potential of being widely applicable, in particular in industrial contexts 
where geometries are commonly defined by NURBS. Simulations indicate that IGS is comparable 
to other widely used methods, such as TPS, in cases where the latter is applicable. Moreover, IGS 
avoids the use of complex meshing procedures since the geometry of the surface is directly used as 
a data of the problem. IGS is also computationally efficient; this is due to the fact that the NURBS 
basis functions have compact support and thus the matrices involved in the computations are usu¬ 
ally very sparse. These are very convenient features, especially if the number of basis functions is 
relatively large (e.g. in the case of a complex surface) and/or when there are many data points. 

A differential penalization is equivalent to the assumption that the function to be estimated 
must be close to the kernel of the penalization operator (Green and Silverman, 1993). In the case 
where the observed physical phenomenon is driven by a known PDE, it would be convenient to 
use a penalization operator related to the PDE (Azzimonti et ah, 2015). Since IGS is based on 
IGA, it would be easy to implement other kind of penalization. IGS also allows the application of 
different kind of boundary conditions, which can be useful in practical applications (Wood et ah, 
2008; Sangalli et ah, 2013; Azzimonti et ah, 2014). In addition, IGS can straightforwardly deal with 
data in one, two, or even three dimensions. Indeed, all the results presented here can be extended 
to any lower dimensional manifold defined by NURBS. Moreover, as shown in the application, 
functionals of an estimated field, as the aerodynamic force, can be easily computed. In addition, 
IGS is very flexible and allows local refinement of the basis functions, which can be needed when 
the distribution of the data points is not uniform. Einally, we remember that it is also possible to 
extend this model to take account of spatially varying covariates on the manifold, using a generalized 
additive framework. 


Acknowledgments 

M. W. and L. D. contributed equally to this paper. The authors aknowledge Prof. Alfio Quarteroni 
and Prof. Eabio Nobile (EPP Lausanne), Denis Devaud (ETH Zurich), and Lionel Wilhelm (EPE 
Lausanne) for very insightful discussions and suggestions. M. W. would like to thank Prof. Yves 
Tille (Universite de Neuchatel) for his support. 


Bibliography 

References 

Alfeld, P., Neamtu, M., Schumaker, L. L., 1996. Pitting scattered data on sphere-like surfaces using 
spherical splines. Journal of Gomputational and Applied Mathematics 73 (1-2), 5-43. 

Anderson, J., 2010. Eundamentals of Aerodynamics. McGraw-Hill Education, New York. 

Azzimonti, L., Nobile, E., Sangalli, L. M., Secchi, P., 2014. Mixed finite elements for spatial regres¬ 
sion with PDE penalization. SIAM/ASA Journal on Uncertainty Quantification 2 (1), 305-335. 


23 


Azzimonti, L., Sangalli, L. M., Secchi, P., Domanin, M., Nobile, F., 2015. Blood flow velocity field 
estimation via spatial regression with PDF penalization. Journal of the American Statistical 
Association 110 (511), 1057-1071. 

Baramidze, V., Lai, M. J., Shum, C. K., 2006. Spherical splines for data interpolation and fitting. 
SIAM Journal on Scientific Computing 28 (1), 241-259. 

Bartezzaghi, A., Dede, L., Quarteroni, A., 2015. Isogeometric analysis of high order partial differ¬ 
ential equations on surfaces. Computer Methods in Applied Mechanics and Engineering 295, 446 
- 469. 

Beaubier, B., Dufour, J.-E., Hild, E., Roux, S., Lavernhe, S., Lavernhe-Taillard, K., 2014. CAD- 
based calibration and shape measurement with stereoDIC. Experimental Mechanics 54 (3), 329- 
341. 

Brezis, H., 1999. Anayse Eonctionnelle: Theorie et Applications. Dunod, Paris. 

Buja, A., Hastie, T. J., Tibshirani, R. J., 1989. Linear smoothers and additive models. The Annals 
of Statistics 17 (2), 453-510. 

Cottrell, J. A., Hughes, T. J. R., Bazilevs, Y., 2009. Isogeometric Analysis: Toward Integration of 
CAD and FEA. Wiley, Hoboken. 

Craven, P., Wahba, G., 1978. Smoothing noisy data with spline functions. Numerische Mathematik 
31 (4), 377-403. 

Dassi, F., Ettinger, B., Perotto, S., Sangalli, L. M., 2015. A mesh simplification strategy for a 
spatial regression analysis over the cortical surface of the brain. Applied Numerical Mathematics 
90, 111-131. 

Dede, L., Quarteroni, A., 2015. Isogeometric Analysis for second order Partial Differential Equations 
on surfaces. Computer Methods in Applied Mechanics and Engineering 284, 807- 834. 

Duchamp, T., Stuetzle, W., 2003. Spline smoothing on surfaces. Journal of Computational and 
Graphical Statistics 12 (2), 354-381. 

Duchon, J., 1977. Splines minimizing rotation-invariant semi-norms in Sobolev spaces. In: Schempp, 
W., Zeller, K. (Eds.), Constructive Theory of Functions of Several Variables. Springer-Verlag, 
Berlin and Heidelberg. 

Dufour, J.-E., Hild, F., Roux, S., 2015. Shape, displacement and mechanical properties from iso¬ 
geometric multiview stereocorrelation. The Journal of Strain Analysis for Engineering Design 
50 (7), 470-487. 

Ettinger, B., Perotto, S., Sangalli, L. M., 2015. Spatial regression models over two-dimensional 
manifolds. Biometrika. 

Green, P. J., Silverman, B. W., 1993. Nonparametric Regression and Generalized Linear Models: 
A roughness penalty approach. Chapman & Hall/CRC, Boca-Raton. 

Hastie, T. J., Tibshirani, R. J., 1990. Generalized Additive Models. Chapman & Hall, CRC Press, 
Boca-Raton. 


24 



Hughes, T. J. R., Cottrell, J. A., Bazilevs, Y., 2005. Isogeometric Analysis: CAD, Finite Elements, 
NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and 
Engineering 194 (39-41), 4135 - 4195. 

Lai, M. J., Schumaker, L. L., 2007. Spline Functions on Triangulations. Cambridge University Press, 
Cambridge. 

Marra, G., Wood, S. N., 2012. Coverage properties of confidence intervals for generalized additive 
models components. Scandinavian Journal of Statistics 39, 53—74. 

Nocedal, J., Wright, S. J., 1999. Numerical Optimization. Springer-Verlag. New York. 

Piegl, L., Tiller, W., 1997. The NURBS Book. Springer-Verlag, New York. 

Quarteroni, A., 2015. Numerical Models for Differential Problems. Springer-Verlag, Milan. 

Quarteroni, A., Valli, A., 1994. Numerical Approximation of Partial Differential Equations. 
Springer-Ver lag, Heidelberg. 

Ramsay, J. O., Silverman, B. W., 2005. Functional Data Analysis. Springer-Verlag, New York. 

Ramsay, T. O., 2002. Spline smoothing over difficult regions. Journal of the Royal Statistical Society: 
Series B (Statistical Methodology) 64 (2), 307-319. 

Sangalli, L. M., Ramsay, J. O., Ramsay, T. O., 2013. Spatial spline regression models. Journal of 
the Royal Statistical Society: Series B (Statistical Methodology) 75, 681-703. 

Stoker, J. J., 1989. Differential Geometry. Wiley, Hoboken. 

Tagliabue, A., Dede, L., Quarteroni, A., 2014. Isogeometric Analysis and error estimates for high 
order Partial Differential Equations in fluid dynamics. Computers & Eluids 102, 277-303. 

Wahba, G., 1981. Spline interpolation and smoothing on the sphere. SIAM Journal on Scientific 
and Statistical Computing 2 (1), 5-16. 

Wahba, G., 1990. Spline Models for Observational Data. Society for Industrial and Applied Math¬ 
ematics, Philadelphia. 

Wilhelm, M., 2013. Generalized spatial regression with differential penalization. Master’s thesis, 
Ecole Polytechnique Federate de Lausanne. 

URL //infoscience.epf1.ch/record/188219 

Wood, S. N., 2006. Generalized additive models: an introduction with application in R. Chapman 
Sz Hall/CRC, Boca-Raton. 

Wood, S. N., 2011. Fast stable restricted maximum likelihood and marginal likelihood estimation 
of semiparametric generalized linear models. Journal of the Royal Statistical Society: Series B 
(Statistical Methodology) 73 (1), 3-36. 

Wood, S. N., 2015. mgcv: Mixed GAM Computation Vehicle with GCV/AIC/REML smoothness 
estimation. R package version on 1.8-6. 

Wood, S. N., Bravington, M. V., Hedley, S. L., 2008. Soap film smoothing. Journal of the Royal 
Statistical Society: Series B (Statistical Methodology) 70 (5), 931-955. 


25 


